function osl_test_script_oil(testdatadir,testoutputdir,S)

global OSLDIR;

forward_meg=S.forward_meg;
work_in_pca_subspace = S.work_in_pca_subspace;
if S.do_hmm ==1; error('hmm not supported by oil.');end

testoutputdir_full=[testoutputdir '_script_oil_rest' '_' regexprep(S.forward_meg, ' ', '') '_' S.recon_method];
%runcmd(['rm -rf ' testoutputdir_full]);
mkdir(testoutputdir_full);

testplotsdir=[testoutputdir_full '/plots'];
mkdir(testplotsdir);

%%

testdatadir=[testdatadir '/oil_test_data/'];

spm_files = {
    
[testdatadir 'tAdspm8_ctf_3-state_ec_1'];
[testdatadir 'tAdspm8_ctf_3-state_eo_1'];
[testdatadir 'tAdspm8_ctf_3-state_mov_1'];
[testdatadir 'tAdspm8_ctf_3-state_ec_2'];
[testdatadir 'tAdspm8_ctf_3-state_eo_2'];
[testdatadir 'tAdspm8_ctf_3-state_mov_2'];
[testdatadir 'tAdspm8_ctf_3-state_ec_3'];
[testdatadir 'tAdspm8_ctf_3-state_eo_3'];
[testdatadir 'tAdspm8_ctf_3-state_mov_3'];
[testdatadir 'tAdspm8_ctf_3-state_ec_6'];
[testdatadir 'tAdspm8_ctf_3-state_eo_6'];
[testdatadir 'tAdspm8_ctf_3-state_mov_6'];
[testdatadir 'tAdspm8_ctf_3-state_ec_7'];
[testdatadir 'tAdspm8_ctf_3-state_eo_7'];
[testdatadir 'tAdspm8_ctf_3-state_mov_7'];
[testdatadir 'tAdspm8_ctf_3-state_ec_8'];
[testdatadir 'tAdspm8_ctf_3-state_eo_8'];
[testdatadir 'tAdspm8_ctf_3-state_mov_8'];
[testdatadir 'tAdspm8_ctf_3-state_ec_9'];
[testdatadir 'tAdspm8_ctf_3-state_eo_9'];
[testdatadir 'tAdspm8_ctf_3-state_mov_9'];
[testdatadir 'tAdspm8_ctf_3-state_ec_10'];
[testdatadir 'tAdspm8_ctf_3-state_eo_10'];
[testdatadir 'tAdspm8_ctf_3-state_mov_10'];

};


%% Beta Band ICA - like Brookes PNAS 2012 

load([OSLDIR '/savedState']);
defaultStream = RandStream.getDefaultStream;
defaultStream.State = savedState;
initGuess = orth(randn(3559,25));

oil = [];
oil.paradigm = 'rest';
oil.source_recon = [];
oil.source_recon.D_continuous = spm_files(2:3:end);                                      
oil.source_recon.mri = {'';'';'';'';'';'';'';'';};
oil.source_recon.conditions = {'Undefined'};
oil.source_recon.time_range='';                                       
oil.source_recon.modalities = {'MEGGRAD'};
oil.source_recon.fid_label.nasion='nas';
oil.source_recon.fid_label.lpa='lpa';
oil.source_recon.fid_label.rpa='rpa';
oil.source_recon.gridstep = 8; % normally 6 mm
oil.source_recon.dirname = [testoutputdir_full '/oil_example_13-30Hz_ctf_eyes_open_data_8subs'];
oil.source_recon.freq_range=[13 30];                                       
oil.source_recon.forward_meg =forward_meg;
oil.source_recon.pca_dim = 0;
oil.source_recon.work_in_pca_subspace = work_in_pca_subspace;
oil.source_recon.force_pca_dim = 0;
oil.source_recon.method=S.recon_method;

oil.enveloping.gridstep = 8;
oil.enveloping.window_length = 1;
oil.ica.num_ics = 25;
oil.ica.temp_or_spat = 'temporal';
oil.ica.initGuess = initGuess;
oil.to_do = [1 1 1 1 0 0];
oil.do_plots=0;

oil = osl_check_oil(oil);

oil = osl_run_oil(oil);

vol = 1:5;   osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EO_beta_band_ica_comps_1-5'],'Comps 1-5',[96 99.9]);
vol = 6:10;  osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EO_beta_band_ica_comps_6-10'],'Comps 6-10',[96 99.9]);
vol = 11:15; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EO_beta_band_ica_comps_11-15'],'Comps 11-15',[96 99.9]);
vol = 16:20; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EO_beta_band_ica_comps_16-20'],'Comps 16-20',[96 99.9]);
vol = 21:25; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EO_beta_band_ica_comps_21-25'],'Comps 21-25',[96 99.9]);

%% Plot the Variance Maps for the group of Eyes Open Subjects

Nvoxels = numel(nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(1)) '_recon_noise_corrected_variance_map'],oil.source_recon.gridstep));
varmaps = zeros(Nvoxels,9); c=1;
for i=1:8
    singsub_varmaps_name = [oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(i)) '_recon_noise_corrected_variance_map'];
    var_maps_singsub = nii_quickread(singsub_varmaps_name,oil.source_recon.gridstep);
    varmaps(:,c) = var_maps_singsub/std(var_maps_singsub);
    c=c+1;
end
varmaps(:,c) = mean(varmaps(:,1:c-1),2);

fname=[oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/average_of_sessions_recon_noise_corrected_variance_map_stacked'];
map_name=nii_quicksave(varmaps,fname,oil.source_recon.gridstep,2);

vol = 1:9;   osl_brain_plot(map_name ,'zslice_plot',42,vol,2,'brain', [testplotsdir '/EO_beta_band_variance_maps'],'Eyes Open Variance Maps',[1 99]);


%% Alpha Band Dual Regression Test

load([OSLDIR '/savedState']);
defaultStream = RandStream.getDefaultStream;
defaultStream.State = savedState;
initGuess = orth(randn(3559,25));

oil = [];
oil.paradigm = 'rest';
oil.source_recon = [];
oil.source_recon.D_continuous = spm_files(setdiff(1:24,2:3:23));                                      
oil.source_recon.mri = {'';'';'';'';'';'';'';'';'';'';'';'';'';'';'';'';};
oil.source_recon.conditions = {'Undefined'};
oil.source_recon.time_range='';                                       
oil.source_recon.modalities = {'MEGGRAD'};
oil.source_recon.fid_label.nasion='nas';
oil.source_recon.fid_label.lpa='lpa';
oil.source_recon.fid_label.rpa='rpa';
oil.source_recon.gridstep = 8; % normally 6mm
oil.source_recon.dirname = [testoutputdir_full '/oil_example_8-13Hz_ctf_eyes_closed_vs_movie_data_8subs'];
oil.source_recon.freq_range=[8 13];                                       
oil.source_recon.forward_meg =forward_meg;
oil.source_recon.pca_dim = 0;
oil.source_recon.work_in_pca_subspace = work_in_pca_subspace;
oil.source_recon.force_pca_dim = 0;
oil.enveloping.gridstep = 8;
oil.enveloping.window_length = 1;
oil.ica.num_ics = 25;
oil.ica.temp_or_spat = 'temporal';
oil.ica.initGuess = initGuess;


oil.ica_group_level.group_design_matrix = [
    1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1;...
    1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
    0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ;...
    0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 ;...
    0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 ;...
    0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0;...
    0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 ;...
    0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0;...
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1;...
    ]';

oil.ica_group_level.group_contrast = {};
oil.ica_group_level.group_contrast{1}=[1 0 0 0 0 0 0 0 0 ]';

oil.ica_group_level.use_randomise = 1; 
oil.ica_group_level.demean_designmatrix = 0;
oil.ica_group_level.group_varcope_spatial_smooth_fwhm = 16;
oil.ica_group_level.Npermutations = 2000; % normally 5000
oil.ica_group_level.cluster_threshold=4;
oil.ica_group_level.comps2use=1:oil.ica.num_ics;

oil.do_plots=0;

oil.to_do = [1 1 1 1 1 1];
oil = osl_check_oil(oil);    

oil = osl_run_oil(oil);

vol = 1:5;   osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_ica_comps_1-5'],'Comps 1-5',[96 99.9]);
vol = 6:10;  osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_ica_comps_6-10'],'Comps 6-10',[96 99.9]);
vol = 11:15; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_ica_comps_11-15'],'Comps 11-15',[96 99.9]);
vol = 16:20; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_ica_comps_16-20'],'Comps 16-20',[96 99.9]);
vol = 21:25; osl_brain_plot(oil.ica.results.maps ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_ica_comps_21-25'],'Comps 21-25',[96 99.9]);

vol = 1:5;   osl_brain_plot(oil.ica_group_level.results.tstats_names{1} ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_tstats_comps_1-5'],'Comps 1-5',[96 99.9]);
vol = 6:10;  osl_brain_plot(oil.ica_group_level.results.tstats_names{1} ,'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_tstats_comps_6-10'],'Comps 6-10',[96 99.9]);
vol = 11:15; osl_brain_plot(oil.ica_group_level.results.tstats_names{1},'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_tstats_comps_11-15'],'Comps 11-15',[96 99.9]);
vol = 16:20; osl_brain_plot(oil.ica_group_level.results.tstats_names{1},'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_tstats_comps_16-20'],'Comps 16-20',[96 99.9]);
vol = 21:25; osl_brain_plot(oil.ica_group_level.results.tstats_names{1},'all_views','max_vox',vol,2,'brain', [testplotsdir '/EC_MOV_alpha_band_tstats_comps_21-25'],'Comps 21-25',[96 99.9]);

%% Plot the Variance Maps for the group of Eyes Open Subjects

Nvoxels = numel(nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(1)) '_recon_noise_corrected_variance_map'],oil.source_recon.gridstep));
varmaps = zeros(Nvoxels,9); c=1;
for i=1:2:15
    singsub_varmaps_name = [oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(i)) '_recon_noise_corrected_variance_map'];
    var_maps_singsub = nii_quickread(singsub_varmaps_name,oil.source_recon.gridstep);
    varmaps(:,c) = var_maps_singsub/std(var_maps_singsub);
    c=c+1;
end
varmaps(:,c) = mean(varmaps(:,1:c-1),2);

fname=[oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/average_of_sessions_recon_noise_corrected_variance_map_stacked'];
map_name=nii_quicksave(varmaps,fname,oil.source_recon.gridstep,2);

vol = 1:9;   osl_brain_plot(map_name ,'zslice_plot',42,vol,2,'brain', [testplotsdir '/EO_beta_band_variance_maps'],'Eyes Closed Variance Maps',[1 99]);

%% Plot the Variance Maps for the group of Eyes Open Subjects

Nvoxels = numel(nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(1)) '_recon_noise_corrected_variance_map'],oil.source_recon.gridstep));
varmaps = zeros(Nvoxels,9); c=1;
for i=2:2:16
    singsub_varmaps_name = [oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/session' num2str(oil.concat_subs.sessions_to_do(i)) '_recon_noise_corrected_variance_map'];
    var_maps_singsub = nii_quickread(singsub_varmaps_name,oil.source_recon.gridstep);
    varmaps(:,c) = var_maps_singsub/std(var_maps_singsub);
    c=c+1;
end
varmaps(:,c) = mean(varmaps(:,1:c-1),2);

fname=[oil.source_recon.dirname '/' oil.enveloping.name '/variance_maps/average_of_sessions_recon_noise_corrected_variance_map_stacked'];
map_name=nii_quicksave(varmaps,fname,oil.source_recon.gridstep,2);

vol = 1:9;   osl_brain_plot(map_name ,'zslice_plot',42,vol,2,'brain', [testplotsdir '/EO_beta_band_variance_maps'],'Movie Variance Maps',[1 99]);
